Energy landscape of a Lennard-Jones liquid: Statistics of stationary points 
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Molecular dynamics simulations are used to generate an ensemble of saddles of the potential 
energy of a Lennard Jones liquid. Classifying all extrema by their potential energy u and number 
of unstable directions k, a well defined relation k{u) is revealed. The degree of instability of typical 
stationary points vanishes at a threshold potential energy uth, which lies above the energy of the 
lowest glassy minima of the system. The energies of the inherent states, as obtained by the Stillinger- 
Weber method, approach uth at a temperature close to the mode-coupling transition temperature 



The complex physical behaviour of supercooled hquids 
and glasses has always stimulated a description of such 
systems in terms of dynamical evolution upon a very 
complicated potential energy landscape In partic- 
ular, the presence of many inequivalent glassy minima 
of the potential energy gives rise to a rich pattern of 
activated dynamics and has therefore attracted most of 
the attention of the community in the past years. As a 
consequence, the analysis of the structure of minima has 
become the main focus of the energy landscape approach 
to glasses. Based on the pioneering work of Goldstein 
101 and of Stillinger and Weber [|| , a number of authors 
have investigated the statistical properties of local min- 
ima, which are reached by a gradient descent, starting 
from an equilibrium configuration of the supercooled liq- 
uid state at temperature T [||-^. The energies of the 
inherent minima sampled after the quench reveal two 
characteristic temperatures: At low T the system gets 
trapped in a very small number of basins and falls out 
of equilibrium. Although this transition temperature de- 
pends on the cooling rate, it is very close to the dynam- 
ical transition Tc of mode-coupling theory (MOT) [Q. A 
second higher temperature marks the onset of noncxpo- 
nential relaxation, which has been associated with energy 
landscape dominated dynamics 

A method which does not focus purely on the proper- 
ties of minima, is the instantaneous normal mode (INM) 
approach ^ : the spectrum of the eigenvalues of the Hes- 
sian matrix is computed and averaged over all the config- 
urations with the Boltzmann distribution at temperature 
T. The INM analysis focuses on two points: First, it has 
been suggested that the barrier heights and hopping rates 
can be obtained from the INM spectrum ||^. Second, the 
temperature where the fraction /(T) of negative eigen- 
values of the INM spectrum goes to zero is interpreted as 
the point where the number of directions for free diffusion 
in phase space vanishes, thus giving an estimate of T^, be- 
low which activation remains as the only mechanism of 
diffusion The INM analysis has been criticized [ p^ , 
because equilibrium configurations with unstable direc- 
tions are in general not stationary points of the potential 



energy, even if the force vanishes along the unstable di- 
rections ijll]] . 

In contrast to the INM method, which is an intrin- 
sically thermal approach by sampling configurations ac- 
cording to their Boltzmann weight, we focus here on the 
purely geometric properties of the energy landscape. We 
investigate all the stationary points of the potential en- 
ergy, be they minima, or unstable saddles We clas- 
sify them according to their number of unstable direc- 
tions, (index K), their energy and the smallest eigenvalue 
of the Hessian matrix. Thereby we can address the fol- 
lowing questions: Is there a threshold energy for saddles, 
such that for energies below this threshold, it is very 
unlikely to find saddles and the dynamic relaxation of 
the system requires activation ? Is there a signature of 
the threshold energy in the dynamical behaviour of the 
system ? What is the typical energy difference between 
saddles of index K and K + 1 and can this be taken as 
an estimate of the potential energy barriers ? 

The system under consideration is a binary mixture 
of large (L) and small (S) particles with 80% large and 
20% small particles. Small and large particles only dif- 
fer in diameter, but have the same mass. They interact 
via a Lennard-Jones potential of the form Vaf}{rij) — 
4eQ/3[((TQ/3/rjj)i2 - {(Jais/nj)'^], where n denotes the po- 
sition of particle i (i = l,2...iV) and = \fi — fj\. 
All results are given in reduced units, where a^L was 
used as the length unit and cll as the energy unit. 
The other values of e and a were chosen as follows: 
eLS = 1.5, aLS = 0.8,655 = 0.5, ass = 0.88 The 
systems were kept at a fixed density p « 1.2. Periodic 
boundary conditions have been applied and the poten- 
tial has been truncated appropriately according to the 
minimum image rule . We have applied a truncation 
procedure which on the one hand ensures the potential 
energy to be zero at the cut-off r^ut and on the other 
hand provides a continuous first-derivative of Vap{r) at 
Tcut- Throughout this study we present results for sys- 
tems with iV = 60 particles using Tcut = 1.8. A few 
samples with up to = 120 have been simulated to val- 
idate our results. 
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In order to explore the geometric properties of the 
stationary points of the potential energy, we use the 
following method: We start by equilibrating a random 
initial configuration of N particles at a given temper- 
ature T using a standard molecular dynamics simula- 
tion technique. After equilibration the system evolves 
for a time Trun- To locate a saddle of the poten- 
tial energy close to the equilibrated configuration we 
look for the absolute minima of the modulus square of 
the force. In order achieve this we perform a quench 
on a pseudo-potential energy landscape W(x) given by 
W{x) = VU{x) ■ VU{x), whereas the original potential 
energy is defined by U{x) J2i<i<j<N ^afiirij). Notice 
that all absolute minima of W{x) are stationary points 
of U{x), hence every saddle of U{x) has a well defined 
basin of attraction. The local minima of W{x) do not 
correspond to zeros of the real force. These points are 
frequently sampled; they can easily be distinguished from 
the absolute minima and are excluded from the following 
analysis. 

Given a stationary point we consider the number of 
negative eigenvalues of its Hessian matrix, that is the 
index K of the saddle. We compute the index density 
k = K/{3N) and the potential energy density u = U/N 
of the saddle and plot these values in the (u, k) plane. 
In Fig.l we show the results for many saddles, which 
were sampled with the steepest descent procedure on 
the pseudo-potential W{x) for two different values of the 
temperature. This plot clearly suggests that there is an 
underlying curve k{u), which is independent of temper- 
ature and which encodes a purely geometric feature of 
the landscape. By sampling stationary points at differ- 
ent values of T we are exploring different regions of the 
potential energy surface and thus different portions of the 
same geometric curve k(u). 



0.08 



0.04 







O T=0.5 
X T=2.0 



XX 
X ^ X 

X x»oa«xx X 
XX X 



o< 



O 



)0 



-2.5 



=4.5 -4 -3.5 -3 

potential energy density 

FIG. 1. Index density as a function of the potential en- 
ergy density: sampling of stationary points at two different 
temperatures, T — 0.5 and T = 2.0. 

To further support this conjecture we have sampled 



stationary points at various temperatures and have av- 
eraged all data in two different ways: Firstly, we have 
considered all the stationary points with a given index 
density and we have computed their average potential 
energy (this is possible because for a finite system k can 
only assume discrete values) . We call this procedure geo- 
metric average. Secondly, we consider all the stationary 
points obtained at a given temperature T and we com- 
pute their average index k{T) and energy u{T). Eventu- 
ally, we plot k and u parametrically in T on the {u, k) 
plane. We call this the parametric average. If our sam- 
pling of the saddles is a fair exploration of the underlying 
geometric space, then the two averages must coincide. 
This is what happens, as it is shown in Fig. 2. 



0.12 



0.1 



0.08 



0.06 



0.04 



0.02 



O Geometric average 




♦ Parametric average 













-4.5 



-2.5 



-4 -3.5 -3 

energy density 

FIG. 2. Index density as a function of the potential energy 
density. Average over all the data obtained by sampling at 
T € [0.3,2.0]. The full line is a linear fit of the geometric 
average. 

The averaged data as plotted in Fig. 2 reveal a unique 
function k{u). In other words, if we cut the potential en- 
ergy landscape with a plane of constant energy density 
u = uq, the stationary points on this plane (or within 
a narrow shell around this plane) will be dominated by 
saddles with index density k{uo). Furthermore, fc(u) is 
to a very good approximation linear up to an index den- 
sity of 10% negative eigenvalues. This implies that the 
curve extrapolates to zero at a well defined energy, which 
we call the threshold energy, uth, in analogy with spin- 
glasses In Fig. 3 we present a magnification of the 
last four points of k{u), showing that the linear interpo- 
lation of all the data and the linear interpolation of the 
last four points give the same estimate for the threshold 
energy, that is Uth = —4.55. 

The threshold energy marks the border between the 
saddles-dominated portion of the energy landscape and 
the minima-dominated one. An interesting point is that 
Uth is above the energy of the lowest lying minima we 
find. Mo = —4.65, as shown in Fig. 3 (uq is obtained from 
an extensive search for minima of the potential energy 
using the Stillinger- Weber method). This implies the ex- 
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istence of a finite energy density interval where minima 
are overwhelmingly more numerous than saddles. The 
same phenomenon has been observed in mean-field mod- 
els of spin-glasses jlsl exhibiting one step replica symme- 
try breaking (fRSB). 
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FIG. 3. Magnification of Fig. 2. With the arrows are indi- 
cated the energy uq of the lowest minima found in the system 
and the threshold energy uth- 

In Fig. 4 we plot the lowest eigenvalue Ao of the Hessian 
as a function of the energy density. As expected, Ao ^ 
for u —> Uth, implying that most extrema have a small 
number of eigenvalues close to zero. Somewhat surpris- 
ing is the approximately linear dependence of Ao(u) on 
energy, Aq {utt — u), exactly as in IRSB spin-glasses 
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FIG. 4. Lowest eigenvalue of the Hessian as a function of 
the energy density. Full diamonds are the parametric average 
of all the data obtained at the same values of T as in Fig. 2. 

It is important to relate the geometric features of the 
energy landscape to the dynamical behaviour of the su- 
percooled liquid close to the glass transition. The in- 
dex density vanishes at the threshold so that minima are 
the dominant stationary points below uth- We there- 
fore expect to find a signature of uth in the dynamics 



or, more specifically, a link between the threshold energy 
and the onset of activated dynamics upon cooling. To 
that end we have used our MD simulations at temper- 
ature T to compute the energy density UminiT) of 
local minima reached by gradient descent on the poten- 
tial energy surface. These data are compared in Fig. 5 to 
the difference SiT) =< U/N > (T) - 3T/2 of the aver- 
age potential energy density < U/N > (T), as calculated 
in the MD simulation, and the vibrational energy in the 
harmonic approximation. For a harmonic potential S(T) 
is just the energy of the minimum of the well. In Fig. 5 
we see that S{T) - (T) for T < 1.2 0. This is 
the range of temperatures which is dominated by the en- 
ergy landscape and the timescales for the two processes 
of relaxation - vibrations inside a minimum and hopping 
between different minima - start to separate. Close to 
the glass transition the system falls out of equilibrium, 
as indicated in Fig. 5 by the saturation of both quantities, 
Umin{T) and 5{T). The temperature where this happens 
is known to depend on the equilibration time of the MD 
simulation (see e.g. [||). 




FIG. 5. Triangles represent the average energy of lo- 
cal minima Umin as a function of the temperature of 
the initial MD trajectory. Circles represent the quantity 
5{T) = (U/N){T) - 3/2T. 

Extrapolating that part of the curve corresponding 
to temperatures where equilibration is ensured (see Fig. 
5), one observes that both Umin{T) and S{T) reach the 
threshold energy at a temperature approximately equal 
to the MCT transition ~ 0.44 The crossover from 
a non-activated dynamics above Tc to an activated one 
below Tc , may thus be interpreted in this context as 
a geometric transition from saddle dominated, to minima 
dominated regions of the potential energy landscape. 

We can obtain an estimate of the energy barrier AU = 
l/{3k') 5.0 from the slope k'{u) of the linear function 
p8[ . This estimate has the right order of magnitude, see 
e.g. Note however that we do not know which saddles 
are accessible from a particular minimum via a dynamic 
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trajectory. Furthermore we have not yet computed the 
entropic contribution to the transition rate, which is nec- 
essary for an estimate of the free energy barrier and hence 
for a comparison with the experimental time scales. 

Finally in Fig. 6 , we compare our result for k{u) to 
an INM analysis. For the latter we have computed the 
average fraction of negative eigenvalues for equilibrium 
configurations as a function of their energy (parametric 
average). We see that for all temperatures the INM in- 
dex is higher than the index of genuine saddles at the 
corresponding energy. Hence the average curvature is 
overestimated by INM, as expected pO|. 
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FIG. 6. INM vs saddle index density as functions of the 
potential energy density. T G [0.3,2.0]. 

In this Letter we have presented a new method to char- 
acterize in a quantitative way the geometric structure of 
the potential energy landscape, which goes beyond the 
analysis of local minima. We have shown that the typ- 
ical saddles index k depends linearly on the potential 
energy density and we have defined a threshold energy, 
Uth, where the degree of instability vanishes. Further- 
more, the threshold energy has been related to the MCT 
transition temperature Tc- 

The method of steepest descent on the pseudo- 
potential W{x) is powerful enough to provide detailed 
information on the properties of saddles, so that one 
can estimate not only energy barriers but also the en- 
tropic contribution to the transition rates. Given that 
we know all eigenvalues of the Hessian at the stationary 
points, we can evaluate the partition function at the sad- 
dle, which is needed as an input for the transition state 
theory. So far we have mainly looked at one system size, 
i.e. 60 particles, and have only checked a few results for 
larger systems. It would be interesting to study the ef- 
fects of varying system size on barrier energies, entropies 
and number of unstable directions. Another important 
point is the mutual accessibility of minima and saddles. 
This latter problem may be tractable by a combination of 
both, steepest descent on the potential energy and steep- 



est descent on the pseudo-potential surface. Work along 
these lines is in progress. 
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